Stable three-dimensional spatially modulated vortex solitons in Bose-Einstein 

condensates 
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We present exact numerical solutions in the form of spatially localized three-dimensional (3D) 
nonrotating and rotating (azimuthon) multipole solitons in the Bose-Einstein condensate (BEC) 
confined by a parabolic trap. We numerically show that the 3D azimuthon solutions exist as a 
continuous family parametrized by the angular velocity (or, equivalently, the modulational depth). 
By a linear stability analysis we show that 3D azimuthons with a sufficiently large phase mod- 
ulational depth can be stable. The results are confirmed by direct numerical simulations of the 
Gross-Pitaevskii equation. 
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Recently, a novel class of two-dimensional (2D) spa- 
tially localized vortices with a spatially modulated phase, 
the so called azimuthons, was introduced in Ref. [ij. Az- 
imuthons represent intermediate states between the ra- 
dially symmetric vortices and nonrotating multipole soli- 
tons. In contrast to the linear vortex phase, the phase 
of the azimuthon is a staircaselike nonlinear function of 
the polar angle. Various kinds of azimuthons have been 
shown to be stable in media with a nonlocal nonlinear re- 
sponse [1, 0, 0, [1, 01 . Note, that azimuthon solutions were 
found in Refs. [J, H, Q by using an approximate (though 
rather accurate) variational approach and separation of 
variables. The first example of exact numerically found 
azimuthon solutions was presented in Ref. [7]. 

The aim of this Letter is to present exact numerical 
solutions in the form of spatially localized nonrotating 
and rotating (azimuthon) multipole solitons in the three- 
dimensional (3D) BEC confined by a parabolic trap. We 
numerically show that the azimuthon solutions exist as 
a continuous family parametrized by the angular veloc- 
ity or, equivalently, by the parameter which determines 
the phase modulational depth. To our knowledge, this is 
the first time that solutions with such a nontrivial (3D 
azimuthon) topology can be found with very high (up 
to machine precision) accuracy. Moreover, by means of 
a linear stability analysis, we investigate the stability of 
these structures and show that rotating 3D dipole soli- 
tons (azimuthons with two intensity peaks) are stable 
provided that the number of atoms is small enough and 
the phase modulational depth is large enough. Thus, we 
present the first example of stable rotating dipole soli- 
tons in media with local nonlinearity. The presence of 
the trapping potential plays a crucial role in stabilizing 
the solitons. The results were confirmed by direct nu- 
merical simulations of the 3D Gross-Pitaevskii equation 
(CPE). 

We consider a condensate at zero temperature confined 
in an axisymmetric harmonic trap. The dynamics of the 
condensate is described by the normalized GPE for the 



wave function 

i^ = -A^/;+(a;2 + y2 + j^V)V-a|^/;|V, (1) 

where appropriate dimensionless units are used, <j = ±1, 
where the +(— ) sign corresponds to attractive (repulsive) 
contact interaction. 

Equation ([T]) conserves the 3D norm (the normalized 
number of particles) N = J |?/^p(ir, z-component of the 
angular momentum = Im J [?/^*(r x VxV^)]^(ir, and 
energy 



E = 



j {|VV^P + (a 



2^2 



(7 



(2) 

We will seek solutions of Eq. ^ which are station- 
ary in the frame rotating with the angular velocity uj. 
In cylindrical coordinates (r, z,(/9), such solutions of the 
form 

?/^(r, (/p, t) = $(r, z^(p — uot) exp(— z/it), (3) 

where ji is the chemical potential in the rotating frame, 
satisfy the equation 



r dr dO'^ ^ dO 



(4) 



where = (p — cot, and resolve the variational problem 
SS = for the functional S = E - j^N - coM^. The 
chemical potential in the laboratory frame /ni is related 
to the chemical potential in the rotating frame /i by 



fii = fi- uoM^/N. 



(5) 



In what follows, we restrict ourselves to the case of a 
spherically symmetric trap ((^ = 1) and attractive in- 
teraction (<j = 1). Two-dimensional case corresponding 

^ 1 ("pan-cake configuration") was considered in Ref. 
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FIG. 1: Dependence of the angular velocity cj of a dipole 3D 
azimuthon on the modulational depth p for three values of /i 
and cr = 1 (attractive interaction). 




FIG. 2: (color online) Numerically found stationary localized 
(a), (c) nonrotating (a; = 0) and (b), (d) rotating multipole 
3D solutions of Eq. JB]) with cr = 1 (attractive interaction), 
Q = 1, and /i = 4.5: (a) Dipole; (b) azimuthon with two 
intensity peaks and a; = 0.13, p = 0.43; (c) quadrupole; (d) 
azimuthon with four intensity peaks and u = 0.3, p = 0.22. 
Shown are the isosurfaces of the soliton amplitude. 



0. Equation (|4]) was rewritten in Cartesian coordinates 



-iu ( X- y— I ^ + cr|<l>p<l> 0. 

\ oy ox J 



(6) 



and solved numerically with zero boundary conditions. 
The numerical method that we use to find solutions of 
Eq. (0) was first introduced by Petviashvili [Si (see also 
Ref. ^) and more recently in Refs. [13, El, H H 
14| . Petviashvili employed a stabilizing factor to suppress 
divergence of the iteration procedure (or convergence to 



zero solution). Petviashvili and the authors of Refs. [13, 
EH, El, considered periodic boundary conditions and 
worked in the spectral Fourier representation. We work 
in physical space and use zero boundary conditions so 
that the method should be modified. We describe it in 
the following way. Equation (j6|) can be written as L$ = 
A^($), where stands for the linear part of the equation 
and N{^) is the nonlinear term. At each nth stage of the 
iteration we solve the corresponding linear problem and 
the iteration procedure is 



where the stabilizing factor 

($;,i-iiV($„))' 



(7) 



(8) 



where (u^v) = J uvdxdy, and for the cubic nonlinearity 
A = 3/2. Note, that one can also take (sometimes it 
yields faster convergence) 



s = 



(L-iiV($*),L-iiV($„))^ 



(9) 



with A = 3/4 in Eq. ([7j). It can be seen that the right- 
hand side of Eq. ([7j) has homogeneity zero with respect 
to ^, and this, indeed, prevents the aforementioned di- 
vergence. The convergence can be monitored using the 
value \s — 1|, which should approach zero. 
As an initial guess we use 



(cos mO -\-ipo sin mO) , 



(10) 

on a Cartesian grid, where m is an integer. The param- 
eter po determines the modulation depth of the soliton 
intensity and without loss of generality we can assume 
^ \po\ ^ 1- The ansatz Eq. (p!Q|) describes topology of 
the 3D azimuthon solutions. Note that the case po = 
corresponds to nonrotating multipole 3D solitons (e. g. 
m = 1 to a dipole, m = 2 to a quadrupole etc.), while 
the opposite case \po\ = 1 corresponds to the radially 
symmetric 3D vortices with the topological charge m. 
A detailed study, including a linear stability analysis, of 
3D vortices in BEC with attractive interactions and a 
parabolic trap was performed in Refs. [H,[l7|].The inter- 
mediate case < |po I < 1 corresponds to the rotating 3D 
azimuthons. 

Note, that Eq. (|6]) with a = 1 has localized solutions 
provided that jj^i < 4 + (7 , where jj^i is the chemical 
potential in the laboratory frame. The fundamental soli- 
tons of Eq. ([6]) with <j = 1 and = 1 exist only if /i/ < 3. 
Moreover, the radially symmetric 3D vortices with m = 1 
are stable if 3.72 < /i/ < 5 (for Q = 1) so that in 

the following we consider only the region 3.8 < /i < 5. 

The rate of convergence does not depend on the choice 
of A, a, and po. In all runs we used A = 1^ a = b = 1^ and 
Po between and 1. During the iterations, the parameter 
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FIG. 3: The growth rates 7 as functions of the phase modu- 
lational depth p for different /x. 

p (modulational depth), which is similar to the parameter 
Po in Eq. (p!Q|) , can be introduced (after ehminating the 
constant phase factor in <l>^) in the following way: 

= max |Im^n|/ max |Re (H) 

We start with the case m = 1. The iteration procedure 
monotonically converges to a solution that consists of two 
dipole-shaped structures in the real and imaginary parts 
of ^ (dipole azimuthon, or, azimuthon with two intensity 
peaks). If the amplitudes of these parts are equal (i. e. 
p = 1) we get the radially symmetric 3D vortex with the 
topological charge m = 1. If these amplitudes are differ- 
ent (0 < p < 1), we have a rotating azimuthon; if one of 
the amplitudes is zero (p = 0), we have the nonrotating 
3D dipole soliton. The final value of p depends only (for 
the fixed ja) on the rotational velocity co and does not 
depend on the initial guess of p (for example, the initial 
condition with cj = and p = 1 converges to the non- 
rotating dipole with p = 0). Typically we used a spatial 
grid with a resolution 100^. 

The progressive iterations in the relaxation method 
were terminated when the norm of the residual fell be- 
low 10"'' — 10~^^ (with the corresponding |5 — 1| < 
10~^ — 10~^^). Under this, the relative differences be- 
tween the successive amplitudes lA^+i — An\/An^ where 
An = max \ ^n\^ and modulational depths \pn-\-i — Pn\/Pn 
were less than 10 ~^ and 10~^ respectively. Typically, 
from 9 (for solutions with small p) to 160 (for solu- 
tions with p close to 1) iterations were required for con- 
vergence. Increasing the number of iterations and grid 
points, we were able to find the solutions with accuracy 
up to machine precision. 

By changing lj one can find the azimuthon solution 
with arbitrary p and, thus, there exists a continuous fam- 
ily of azimuthons parametrized by the parameter p. The 
dependence of the angular velocity co of the azimuthon 



with two intensity peaks on the modulational depth p for 
several different ja and a = 1 (attractive interaction) is 
presented in Fig.[Tl In Figs. [2ja) and (b) we demonstrate 
two numerically found examples of the azimuthons (non- 
rotating dipole and rotating azimuthon with two inten- 
sity peaks) for the model described by Eq. ([T]). Similar 
solutions can be found for a = —1 (repulsive interaction) 
and fii > 5. 

The vortex solitons {p = 1) have the maximum angu- 
lar momentum M^, while the dipole solitons {p = 0) have 
zero angular momentum. Structures with the intermedi- 
ate values of p carry out the nonzero angular momentum, 
and this leads to soliton rotation. 

Changing topology of the initial guess one can find 
high-order azimuthon solutions. We took Eqs. (p!Q|) with 
m = 2. The iteration procedure began to converge mono- 
tonically to a solution that consists of two quadrupole- 
shaped structures in the real and imaginary parts of <l>. 
If the amplitudes of these parts are equal {p = 1) we get 
the radially symmetric vortex with the topological charge 
m = 2. In the intermediate case (0 ^ p < 1), we have 
the azimuthon with four intensity peaks. The conver- 
gence was controlled by stopping the iteration when the 
value |5 — 1| began to increase. This indicates that the 
iteration procedure jumps off the solution correspond- 
ing high-order azimuthon (and begins to converge to the 
azimuthon with two intensity peaks). Nevertheless, the 
high-order solutions can be found with a rather high ac- 
curacy: typically, one can reach |5 — 1| ~ 10~^. Then, the 
obtained solution can be used as an initial condition in 
the Yang-Lakoba iterative procedure [lH which belongs 
to a family of universally-convergent iterative methods 
and can converge to any nonfundamental (i. e. high- 
order) solution of a given equation provided that the ini- 
tial condition is sufficiently close to that solution. Thus, 
we were able to find high-order azimuthons of Eq. ([1]) 
with very high accuracy. The nonrotating quadrupole 
and rotating azimuthon with four intensity peaks are pre- 
sented in Figs. [Sfc) and (d). 

The physical relevance of any stationary solution, of 
course, depends on whether it is stable. To study the sta- 
bility of the stationary solutions, we represent the wave 
function in the form 

i^{r,t) = [Mr)+e{r,t)]e-'^^*, (12) 

where the stationary solution cpo is perturbed by a small 
perturbation s. In the absence of the radial symmetry, 
the corresponding eigenvalue problem (i. e. after lin- 
earizing and taking £:(r, t) ~ u{r) ex.p{jt)) on a spa- 
tial grid implies a 2N^ x 2N^ complex nonsymmetric 
matrix and, for reasonable N (say, N > 100), represents 
extremely difficult task. Instead, after inserting Eq. (p!2|) 
into Eq. ([T]), we solved the Cauchy problem for the lin- 



4 




FIG. 4: (color online) The top row: unstable evolution of the 
3D nonrotating (u = 0,p = 0) dipole with /x = 4.2. Initial 
state of the dipole is unperturbed. The middle row: unstable 
evolution of the 3D azimuthon with two intensity peaks and 
a; = 0.18, p = 0.3, /x = 4.2. The initial state is perturbed by 
the white noise with the parameter e — 0.005. The bottom 
row: stable dynamics of the 3D azimuthon with two intensity 
peaks and oo = 0.285, p = 0.73, /i = 4.2. The initial state is 
perturbed by the white noise with the parameter e = 0.08. 

earized equation 

de 
ot 

+2|(/;o|'^ + (/^o^* = (13) 

with some initial small perturbations s. The final results 
are not sensitive to the specific form of ^(r, 0). The chem- 
ical potential in the laboratory frame jui is determined 
from Eq. ([5]). If the dynamics is unstable, the corre- 
sponding solutions ^(r, t), undergoing, generally speak- 
ing, oscillations, grow exponentially in time and an esti- 
mate for the growth rate 7 can be written as 

where P{t) = J l^^dr, t and At are assumed to be large 
enough. 

In Fig. Owe plot the growth rates 7 as functions of the 
modulational depth p for the azimuthons with two inten- 
sity peaks and several different values fi of the chemical 



potential in the rotating frame. The linear stability anal- 
ysis shows that for solutions with p < 0.7, the growth 
rate of perturbations 7 7^ for all ja so that there is 
no stability region. In particular, all nonrotating dipoles 
are unstable. The growth rate 7 decreases as ja increases 
and can be very small if ja is close to 5. The situation, 
however, changes when the modulational depth exceeds 
the critical value Pcr ~ 0.7. In this case, the growth rate 
of perturbations 7 falls to zero, the stability window ap- 
pears and the azimuthons with p ^ 0.7 are stable. The 
critical value Pcr varies very slightly with changing the 
chemical potential ja. All high-order azimuthon solutions 
turn out to be unstable (for the radially symmetric vor- 
tices with p = 1 and the topological charge m > 2 it was 
shown in Ref. 

To verify the results of the linear analysis, we solved 
numerically the dynamical equation ([T]) initialized with 
our computed solutions with added Gaussian noise. The 
initial condition was taken in the form (/?o[l+e^(r)], where 
(po{r) is the numerically calculated solution, ^(r) is the 
white gaussian noise with variance cr^ = 1 and the pa- 
rameter of perturbation e = 0.005 ^0.1. The unstable 
dynamics of the nonrotating dipole {p — 0) with /i = 4.2 
is illustrated in the top row of FigHl In the middle row 
we present unstable evolution of the rotating azimuthon 
with two intensity peaks and \i — 4.2, a; = 0.18, p = 0.3. 
A slight stochastic perturbation with e = 0.005 was ap- 
plied at t = 0. The azimuthon lost its shape after one 
rotational period. Stable evolution of the azimuthon with 
/i = 4.2, = 0.285 and p = 0.73 (i. e. in the region of 
stability) is shown in the bottom row of FiglH The initial 
state of the azimuthon is perturbed by a rather strong 
noise (e = 0.08). The period of rotation of the azimuthon 
is T = 27r/c<; ~ 22. The azimuthon cleans up itself from 
the noise and survives over many dozens of the rotational 
periods. 

In the present paper we restricted ourselves to the case 
of a spherically symmetric trap (^7 = 1). A linear sta- 
bilty of the radially symmetric 3D vortices (i. e. az- 
imuthons with p = 1 in our case) was performed for dif- 
ferent Q in Ref. [17]. By analogy with the results of Ref. 

for the case 1] <C 1 ("cigar-like configuration") one 
could expect that the quasi- ID azimuthon solitons with 
p < 1 are tightly confined in the corresponding cigar- 
like trap, which suppresses their destabilization. On the 
other hand, the case 1 ("pan-cake configuration") 

implies nearly 2D azimuthon solitons, with the largest 
energy, which are most prone to instabilty. 

In conclusion, we have presented exact numerical so- 
lutions in the form of spatially localized 3D nonrotating 
and rotating (azimuthon) multipole solitons in the three- 
dimensional BEG confined by a parabolic trap. We have 
shown that the 3D azimuthon solutions exist as a con- 
tinuous family parametrized by the angular velocity (or, 
equivalently, the phase modulational depth). By a linear 
stability analysis we have shown that 3D azimuthons with 
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a sufficiently large modulational depth can be stable. The 
results are confirmed by direct numerical simulations of 
the Gross-Pitaevskii equation. 

The author thanks A. I. Yakimenko and Yu. A. Zal- 
iznyak for discussions. 
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